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Abstract 

In this article we review recent results concerning numerical simulation 
of shock waves using spectral methods* We discuss shock fitting techniques as 
well as shock capturing techniques with finite difference artificial 
viscosity. We also discuss the notion of the information contained in the 
numerical results obtained by spectral methods and show how this information 
can be recovered. 
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Introduction 


In the last decade spectral methods have been used very successfully in 
the numerical simulations of incompressible flows. Spectral methods have also 
emerged as a major tool in computational meteorology. This has led many 
researchers to look into the possiblity of applying spectral methods to 
simulate compressible flows that are of interest to aeronautical engineers. 
The aim of this article is to give a brief review of the major developments in 
this field in the last few years. In particular we would like to discuss the 
notion of the information that is contained in the numerical result. We argue 
that spectral methods yield more information about the exact solution than low 
order methods. This information is hidden in the form of numerical 
oscillations when the exact solution is discontinuous or contains extreme 
gradients. The structure of these wiggles depends on the nature of the 
discontinuity and, in some cases, a very accurate solution can therefore be 
extracted. 


2. Spectral Methods 

There are basically two steps in obtaining a numerical approximation 
u^(x) to a solution u(x) of a differential equation. First, an appropriate 
finite or discrete respresentation of the solution must be chosen. This may 
take the form of an interpolating function between the values u(xj) at some 
suitable points Xj or a series coefficient in the finite representation 

N 

u N (x) = I \ 4> k (x) 
k=0 


( 2 . 1 ) 
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with given expansion functions <|> k (x). The second step is to obtain equations 
for the discrete values u N (xj) or the coefficients a k from the original 
equations* This second step involves finding an approximation for the 
differential operator in terms of the grid point values of u N or, 
equivalently, the expansion coefficients. For example, the pseudospectral 
Chebyshev approximation to the equation 


u. = u 


u(x,0) = 


u Q (x), 


x | < 1 


u( 1 , t) = h(t) 


(2.2) 


is obtained in the following manner. For a given time t we assume that 
{u^(Xj,t)} is known where x^ = cos . We then interpolate these values 
to get 

N 

(2.3) 


N 

u N (x,t) = l u (x ,t) g. (x) 
j=0 J J 


where 


g d (x> 


(- 1 ) 

N 


j+1 (1 - x 2 )T^(x) 

2 c.(x - x.) ’ 


C 0 = °N = 2 

Cj = 1, 0 < j < N. 


Note that 


SjV ■ V- 


Equivalent ly , s ince 


V 0 


9 N 

= - y 

N L n 
n=0 


T (x.) T (x) 
n j n 

c 

n 


where T n (x) = cos(n cos -1 x) is the Chebyshev polynomial of degree 


n, one 
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u (x,t) - l a T (x) 
N ~ n n 

n=0 


-V 1 *' 3 " O.-.N. 

n C n W j=0 3 C j 


(2.4) 


The next step is to differentiate (2.3) to get the system of ordinary 
differential equations 


3u v.(x, >t) N 

= ^ 0 u M (x r t)g j (x k ): 


j = 1 , • • • ,N 


(2.5) 


3 — <x 0» t) = h ' (t) 


or using (2.4) 


where 


3u vt(x, ,t) N 

^-3 = l a T’ (x, ) 

3t nSo n n ^ 


'TT ( *0' t) ■ h ' (t) 


I b n W- 3 ' l > — "» 
n=0 


( 2 . 6 ) 


3 N “ °» b N-l 2N a N’ b n c b n+2 + 2(n+1)a n+l * 


Equations (2.5) and (2.6) are, in fact, identical. Equation (2.5) points out 
the possibility of applying the pseudospectral Chebyshev method by mulitplying 
the vector u(xj,t) by the matrix gj(x n ) whereas the asymptotically 
efficient implementation of (2.6) is by using a Fast Fourier Transform. 
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In general, consider the system of equations 


u t = L(u) 


u(t=0) = Uq, 


(2.7) 


where L is a nonlinear operator that involves only spatial derivatives. In 

spectral methods we define a finite dimensional subspace Bjj which is the 

space of polynomials (or trigonometric polynomials) of degree N, and a 

projection operator P N that maps the original space to b n . An example of 

such a Pjj is given in (2 • 3) • In fact, given a function f(x), -1 < x < 1, 

N 

then (2.3) defines P f = l f(x.)g.(x). 

j=0 J J 

We then seek a solution u N belonging to B N such that 


9u N 

at “ p n l ^ u n^» 


U N (t=0) = P N V 


( 2 . 8 ) 


For a more complete description of spectral methods we refer the reader to 
[3], [6]. 

Spectral methods are global in nature, i.e., in order to get an expression 

g 

for 3 jc U N we use t * ie points x^, k = 0,*«»,N (see (2.5)). Together 

with the choice of the points xjj this explains their high order accuracy. 
The accuracy of spectral methods depends on the total number of points N, and 
the number of smooth derivatives of u. For smooth flows, great savings of 


computer storage and time is gained by using spectral methods since only a 
small number of grid points is required to get the same accuracy obtained by 
other methods. 


3. Spectral Methods and Shock Haves 

The use of any formal high order method for the numerical simulation of 
flows with shocks poses theoretical and practical problems. The error 
estimates obtained for spectral methods depend on the smoothness of the 
solution and it is not clear at all that any degree of accuracy can be 

achieved for discontinuous solutions. On the one hand, it has been proven 
that for linear problems, high accuracy can be maintained within spectral 
methods far away from the discontinuity; on the other hand, it may be thought 
that for nonlinear problems the overall accuracy in the presence of 
discontinuities is limited to first order. However, in [10] Lax has argued 
that more information about the solution is contained in high resolution 
schemes, even in the nonlinear case. In fact, Lax has shown that the e- 

capacity of the set of approximate solutions is closer to the e-capacity of 

the set that includes the projections of exact solutions if the numerical 

scheme is a high order scheme. Typically, when a spectral method is used to 
simulate flows with shocks it yields an oscillatory solution. The 
oscillations are global, that is they occur not only in the neighborhood of 
the shock but all over the flow field. Several methods of overcoming these 
oscillations were suggested. Historically, the first attempts to get 
nonoscillatory results concentrated on using finite difference type artificial 
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dissipation. Taylor, et al. [15] used the method of Boris, and Book of adding 
diffusion and antidiffusion terms for some model problems* Sakell [12] has 
checked a version of the Von Neumann-Richtmyer artificial dissipation for the 
wedge flow problem. Cornille [2] has used a version of the Lax-Wendroff 
scheme with inherent dissipation. Zang and Hussaini [16] simulated slightly 
viscous flows and treated the viscosity term by finite differences. Two real 
life flows were simulated using the above ideas. Reddy [11] introduced 
Fourier representation in the azimuthal direction in the three-dimensional 
Navier-Stokes code of Pulliam and Steger. In this problem there is enough 
dissipation coming from the discretization in the other directions. Reddy 
reports substantial improvement over the finite difference code. Streett [14] 
simulated transonic flow around an airfoil. His code is a full potential 
algorithm with retarded density. His results indicate that for subsonic 
flows, spectral methods are superior to the finite difference codes, whereas 
for transonic flow they are comparable. The results obtained by these methods 
indicate that a highly structured flow field is well— represented along with 
the front of the shock. However, the shock profiles are smeared and the 
accuracy in the smooth part of the flow is perhaps no longer spectral. 

A different approach advocated first by Hussaini, Salas and Zang [9] is to 
fit the shock. This approach has been used to simulate various physical 
problems, most of them concerned with shock wave interactions. Since they were 
interested in the behavior of the flow on only one side of the shock, a 
coordinate transformation was employed so that the shock wave became a 
coordinate boundary. The Rankine-Hugoniot conditions were used both to 
determine the flow variables immediately upstream of the shock and to 



-7- 


determine the shock position* Since all the physical quantities on the 
downstream of the shock were prescribed the flow variables on the upstream 
side were obtained from the Rankine-Hugoniot relations. Note that the shock 
boundary is supersonic and therefore all the quantities must be specified and 
no special boundary treatment is necessary. The fluid motion was modeled by 
the two-dimensional Euler equation in nonconservation form. Also a spectral 
filtering in which the high modes were filtered every fifty time steps was 
employed to avoid nonlinear instability. Beautiful results were obtained for 
various shock interactions and for the blunt body problem. 

In the third approach proposed in a forthcoming paper by Abarbanel and 
Gottlieb, the oscillations are being used to recover accurate information 
about the solution. Oscillations may arise from different sources; e.g. , 
incorrect treatment of the boundaries in hyperbolic systems; nonlinear 
instabilities, etc. Usually these oscillations build up and finally cause 
explosive instabilities. One interesting class of numerical oscillations 
occur when flows with extreme gradients or local discontinuities are 
simulated. This type of oscillations does not cause instabilities even after 
many time steps. It has been observed (see [7]) that the wiggles are caused 
by the fact that the mesh is not fine enough to resolve the sharp gradients. 
In the case of a finite gradient a local refinement of the mesh often gets rid 
of the wiggles. For a very impressive demonstration of this fact, see [17]. 
Of course for a shock wave, no refinement of the mesh can remove the 
oscillations. 

To better understand the origin of the oscillatory solution, consider the 


model equation 
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u(x,0) = H(x,x^) 

where H(x,x^) is the Heaviside function 

H(x,x^) =0 x < x £ 

H(x»x^) =1 x > x^ 

= cos ( 1 +^/ 2 ), i integer. 


(3.1) 


When (3.1) is discretized by the pseudospectral Chebyshev method we get as the 
initial condition 

N 

u N (x,0) = S(x,x £ ) = l A k T k (x) (3.1a) 

k=0 

where T^(x) is the Chebyshev polynomial of order k, and 
Aq — (^ + V 2 ) > A^ — sin ir (A + V 2 ) 

\ = ¥ sin “If + V2 )/sin 7^- , 1 < k < N-l. 

11 

N 

s C x j » X A ) ■ h(xj,x 4 ), 


At the grid points, 


X j 


cos 
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Thus, no oscillations occur. However, after the numerical solution is 

convected by equation (3.1), it becomes oscillatory. This is because 
initially it is oscillatory between the grid points (see Fig. 1). Observe 
that the oscillations disappear when the discontinuity is exactly in the 
middle between two grid points. This demonstrates the fact that the structure 
of the oscillations provides information about the position and magnitude of 
the shock. 



- 1.0 -.8 -.6 -.4 -.2 0 .2 .4 .6 .8 1.0 

X axis. 


Figure 1 
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In general, consider (2.7) - (2.8) where now L is a linear operator 
and uq is discontinuous. From the last example it is clear that u N does 
not approximate well Pjj. u since Pjj u coincides with u at the grid 
points. We introduce an auxiliary equation 


3v 

8t 


Lv 


v(t=0) 


? N V 


(3.2) 


For fixed N, v is a smooth function in contrast to the solution u of 
(2.7). We argue that u N approximates (at the grid points) v rather than 
u. In fact from (2.8) and (3.2) one gets 


3t ( U N “ P N V ) P N “ P N V ) + P N L ^ P N V ” V ) 

( U N “ P N v)(t=0) = 0. 


(3.3) 


U N " P N V “ f t ex P P N LP N (t - T )][ p N L( p N v(t) - v(t))]dx. 

0 

The operator exp P N LP N (t - t) is bounded. This is, in essence, the notion 
of stability. The term 


P » L ( P N T 


-Il- 


ls small because v is a smooth function. This shows that u^ approximates 
Pjj v, hence at the grid points u ^ approximates v. 

In the last example we have demonstrated the fact the v is, in general, 
oscillatory. It is therefore no surprise that u N is oscillatory. It is 
also clear that the structure of the oscillations may be used to extract a 
better approximation to u. 

We will demonstrate now the possibility of extracting information from an 
oscillatory solution even in the nonlinear case. The physical problem is the 
well-known wedge flow. A plate is inserted in a uniform flow, and an oblique 
shock develops. The time dependent Euler equations in two-space dimensions 
were discretized by the pseudospectral Chebyshev method in space with 
a 9x9 grid and a modified Euler scheme was used for the time discretization 
(see. [5]). Since we are interested in the steady state only, the accuracy of 
the time integration is of no importance. In order to be sure that a steady 
state is reached the code was run until all the physical quantities did not 
change to 11 significant figures over a span of 100 time steps. The values of 
the density in the steady state at the grid points together with the grid 
points themselves are given in Fig. 2. 


p 

Y 


1 .862 

1.851 

1 .869 

1.871 

1.837 

1 .865 

1.892 
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1.878 

1 . 
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1.870 
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1.820 
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1.803 

1.759 
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1.852 
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2.387 

3.133 

3.375 

3.224 

3.054 

3.002 
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3.288 
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3.136 

3.187 
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3.087 
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3.024 

3.013 

3.016 

0 

X 
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.961 
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Figure 2 
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Note that at the stations: xq = 1; x^ a • 96 1 9 ; X£ = *85355, the jump takes 

place between the grid points y = *3086 and y = .5, whereas the 

corresponding correct shock location is y = .434 for x Q , y » .417 for xj 
and y - .370 for X£. Note also that the oscillatory behavior of the density 
is very similar to the behavior of P N v, the solution of (3.2) at the grid 
points (see Fig. 1). 

We therefore fit a step-function of the form dj + d£ S(y,y^) where 
SCy.y^) is defined in (3.1a) to the numerical results p(y) in Fig. 2, at 
any station Xj, regarding dj, d£ and l as unknowns. This yields three 


d l f 0 + d 2 f l " S 1 


d l f l + d 2 f 2 ~ S 2 


(3.4) 


d l f 4 + d 2 f 3 “ S 3 


f 0 = N; f i = l S(y ,y ) -i ; 

j=0 J * c j 


f 2 = I s(y,y) 2 -I; 
j=0 J * C j 


f 3 = k It s ( y j> y J ^ * ^4 j 0 H ( y j. y J ^ ; 


s i = 1 p(yJ — ; 

j-o J C j 


S 2 'jL p(y J )S(y J* y *J ^ ; 


Equation (3.4) yields the following nonlinear equation for the shock location 



(3.5) 


Surprisingly, from (3.5) we recover the correct location of the shock at each 

x-station within the fourth significant digit. In this sense the information 

is indeed hidden in the form of oscillations. 

It should be noted that in (3.4) we do not use the point values of 

p (y) but rather the quantities Sj, S 2 , S 3 which are equivalent to the 

g 

integral of p(y) against 1, S(y,y^) and S(y,y^). If p(y) approxi- 

mates well the first N modes of the solution p £xt (y), then 


1 

/ 


-1 


( p(y) " p ext (y) ) 


-Mi L_ = 0 

/l - y 2 


where <{> (y) 


is either 1 or S(y,y^) 


or ^(y.y*). 


This may be the reason 


for the highly accurate values of the location of the shock obtained by (3.4) . 

Finally, we would like to describe another way of recovering correct point 
values from an oscillatory approximation. For simplicity we consider the 
spectral Legendre method although this idea has been generalized to other 
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spectral methods. Our approach is motivated by the work of Mock and Lax (see 

[ 10 ]). 

Suppose that f(x) is a c“ function at |x| < 1 except for one point 
of discontinuity. Suppose also that f(x) has the following expansion in 
terms of the Legendre polynomials 


00 


f(x) - la. P (x) 
k=0 K 


and that 


N 


k=0 K k 


Even for large N, %(x) is an oscillatory function. Let y be a point such 
that f(x) is C in the interval y-e < x < y+e. Let 




C 2 ) q f (2k+l)P (0)P. (5) 
k=0 k R 


'l'(x) = 


5 < 1 


id > i 


5 


It is clear that 


, 1 1 

J f N (x)i|/(x)dx - J f(x)i(»(x)dx + / (f - f)ij>dx. 

-1 -1 -1 N 


The function ^(x) has the expansion 


+(x) - l b k P (x) 
k=0 
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and since ^(x) has q-1 continuous derivative the function tpjj(x) 


f N (x) - I b. P k (x) 
N k=0 K K 


approximates ^ with high accuracy. Moreover, since is a polynomial 

of degree N 


1 1 

/ (f N - f)tdx - /^(f N - f)0 - ^ N )dx < If - f N « " + N « 
„ .. Ii|^ q-1 ^l 

^ 1C 11 11 , * 


The last estimate can be found in [1]. 


It is therefore clear that 


/ f^ dx «* / f^ dx + E^ 


where is small. Moreover, 


1 1 9 p 

/ f(x)t(x)dx = / f(y+e?)(l-£ ) q £ (2k+l)P. (0)P, (S)d£. 

-1 -1 k=0 


g(5) - f(y + e5)(l- 5 ) . 


g(£) is a C“ function for |£| <1 and therefore has a rapidly converging 
expansion of the form 
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g(5) = l \ P k (0. 
k=0 


Therefore 


IP p 

/ g(0 I (2k+l)P k (0)P k (5)d? = l c k P fe (0) 

-1 k=0 k=0 

00 

= g(0) - l c k = f (y) + E 2 
k=p+l 


This shows that 

/ f N * dx 

approximates f(y) to a high order of accuracy* This filter had been 
successfully used by Gottlieb and Gruberger for several problems* 

In conclusion we have demonstrated that numerical solutions obtained by 
spectral methods contain information about the correct solution that may be 
extracted to yield a high order approximation in the regular sense* 
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